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Abstract. We combine a Monte Carlo radiative transfer code with an SPH code, so that - assuming thermal equilibrium - we 
can calculate dust-temperature fields, spectral energy distributions, and isophotal maps, for the individual time-frames generated 
by an SPH simulation. On large scales, the radiative transfer cells (RT cells) are borrowed from the tree structure built by the 
SPH code, and are chosen so that their size - and hence the resolution of the calculated temperature field - is comparable with 
the resolution of the density field. We refer collectively to these cubic RT cells as the global grid. The code is tested and found 
to treat externally illuminated dust configurations very well. However, when there are embedded discrete sources, i.e. stars, 
these produce very steep local temperature gradients which can only be modelled properly if - in the immediate vicinity of, and 
centred on, each embedded star - we supplement the global grid with a star grid of closely spaced concentric RT cells. 
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1. Introduction 

Smoothed Particle Hydrodynamics (SPH) (Lucy 1977; 
Gingold & Monaghan 1977) is a Lagrangian computational 
method which invokes a large array of particles to describe a 
fluid, assigns properties such as mass, position and velocity to 
each particle, and estimates intensive thermodynamic variables 
like density and pressure (and their derivatives) using local av- 
erages (for reviews see Benz 1990, 1991; Monaghan 1992). 
The main advantage of this method is that, in principle, there 
are no limitations on the geometry of the system or how far it 
can evolve from the initial conditions. 

Because of the large computational cost incurred when ra- 
diative transfer is modelled fully in three dimensions, many 
previous treatments have been limited to one or two dimensions 
(e.g. Efstathiou & Rowan-Robinson 1990, 1991; Men'shchikov 
& Henning 1997, Zucconi et al. 2001), and radiative transfer 
has only recently been included in SPH simulations, where a 
fully three-dimensional treatment is essential (Kessel-Deynet 
& Burkert 2000, Oxley & Woolfson 2003, Whitehouse & Bate 
2004, Viau et al. 2005). The Monte Carlo approach for equilib- 
rium radiative transfer (Lefevre et al. 1982, 1983; Wolf et al. 
1999) is well suited to systems with arbitrary geometries, but, 
as with any Monte Carlo method, it is computationally expen- 
sive, especially when iteration is required. Recently, Bjorkman 
& Wood (2001), extended an idea by Lucy (1999) and pro- 
posed a method to avoid iteration, by remitting photons as soon 
as they are absorbed, with a frequency distribution adjustment 
technique. This method avoids iteration, and it is fast when 
compared with the traditional Monte Carlo radiative transfer 
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methods. It has been applied successfully to a variety of prob- 
lems, such as protoplanetary discs (e.g. Wood et al. 2002), 
prestellar cores (e.g. Stamatellos et al. 2004) and Class I ob- 
jects (e.g. Whitney et al. 2003). 

In this paper, we develop and test a method which uses 
Monte Carlo radiative transfer with frequency distribution ad- 
justment, to calculate the dust-temperature field, the spectral 
energy distribution (SED), and isophotal maps, for individual 
time-frames from SPH simulations. We use the SPH tree to 
construct a grid of cubic radiative transfer cells (hereafter RT 
cells) with which the photons interact; we refer to this grid as 
the global grid. The global grid obtained in this way is simi- 
lar to the tree-structured adaptive grid invoked by Kurosawa & 
Hillier (2001), but because it borrows the tree structure already 
derived as part of the SPH code, it minimizes the computa- 
tional burden that will be added to the hydrodynamic simula- 
tions when the two methods are combined. A similar method 
has been used by Oxley & Woolfson (2003), but they adopt 
an average, wavelength-independent opacity. To account for 
the large temperature gradients that are expected very close to 
stars, we supplement the global grid described above with a 
grid of concentric spherical radiative transfer cells around each 
star; we refer to this supplementary grid as the star grid. 

2. Monte Carlo radiative transfer with frequency 
distribution adjustment 

The Monte Carlo radiative transfer method (e.g. Stamatellos & 
Whitworth 2003) uses a large number of monochromatic lu- 
minosity packets (L-packets) to represent the radiation injected 
into the medium (e.g. stellar radiation, cooling radiation, back- 
ground radiation). Each L-packet is given a random frequency, 
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v, and a random optical depth, r, which determines its free 
path, i.e. how far it travels before interacting with the medium. 
The medium itself is divided into a number of RT cells, each 
with uniform density and temperature. Each time an L-packet is 
scattered, it is given a new random direction, using the Henyey 
& Greenstein (1941) scattering phase function, and a new ran- 
dom optical depth, t. Each time an L-packet is absorbed, it 
raises the temperature of the RT cell in which it is absorbed 
from T to T + AT, and the packet is immediately reemitted 
isotropically, so that the dust is in radiative equilibrium. The 
temperature of the RT cell is computed by equating the total 
absorbed luminosity to the total emitted luminosity. The fre- 
quency of the reemitted L-packet is chosen using a probability 
distribution function (PDF) constructed from the difference be- 
tween the emissivity of the cell before and after absorption of 
the L-packet (Bjorkman & Wood 2001; Baes et al. 2005), 

K v [B v (T + AT)-B v (T)]dv 

p(v)dv = 

J k v [B V (T + AT) - B V (T)] dv 

^ KyB'(T)dV 

J-oo " \ J 

KyB'(T)dV 

where B'(T) = dB v /dT . At the same time the L-packet is given 
a new random optical depth, r. Using the PDF of Eqn. Q, the 
correct temperature distribution and spectrum of the system are 
obtained without iteration, once all the L-packets have been 
propagated through the medium and escaped. The method ac- 
counts for both absorption and scattering, it is robust, it is very 
accurate (see tests by Bjorkman & Wood 2001, Stamatellos & 
Whitworth 2003), it can readily be parallized, and it can treat 
systems with arbitrary geometries. 

2.1. Construction of radiative transfer cells 

The radiative transfer cells must be constructed so that 
their linear size is less than, or on the order of, the lo- 
cal directional temperature- and density- scale-heights (e.g. 
Ajc < MlN{\d£n[p]/dx\-\ \d€n[Tydx\~ 1 }, see Stamatellos 
& Whitworth 2003). This ensures that the temperature and den- 
sity do not vary too greatly between adjacent RT cells. The con- 
struction of the RT cells depends on the specific system under 
study, and, in an hydrodynamic simulation, where the system 
changes from step to step, the RT cells need to be reconstructed 
at every step. Therefore, a robust and efficient algorithm for the 
construction of RT cells is required. 

To construct RT cells we take advantage of the fact that 
our SPH code calculates gravitational forces with the aid of a 
spatial tessellation tree (Barnes & Hut 1986; Hernquist & Katz 
1989). The whole computational domain is contained within a 
cubic cell, the root cell, which is then divided into eight smaller 
cubic cells of equal size. If any of these cells contains more than 
one particle, it is subdivided into eight even smaller cells. This 
procedure continues recursively until the smallest cells each 
contain just one particle or no particles at all. The resulting 
ensemble of cells constitutes the SPH tree. In SPH, the tree is 
used so that distant particles do not contribute individually to 
the local gravitational field, but collectively through the cell in 
which they are contained; the tree is also used to find neigh- 



bouring particles so that smoothed values of thermodynamic 
variables and their derivatives can be evaluated. 

The global grid then comprises the smallest set of RT cells 
which (a) spans the entire computational domain (or at least 
that part of it containing matter) and (b) contains fewer than 
^max SPH particles. If N MAX is chosen to be of the same order 
as N NEm (the number of neighbour particles used for smoothing 
purposes), then the RT cells are necessarily similar in size to 
the resolution of the SPH simulation. Specifically we choose 
N MAX /N NEm ~ 2 to 3, so that each RT cell has linear size ~ 2 h 
(where h is the SPH smoothing length) and contains ~ 0.3 N NElB 
particles. 

This method of constructing RT cells is similar to the 
method of Kurosawa & Hiller (2001) (see also Kurosawa et al. 
2004), but it is implemented using the Barnes & Hut (1986) tree 
structure. Both methods produce RT cells with approximately 
equal mass. It is also similar to the method used by Oxley & 
Woolfson (2003), but they construct RT cells with size 1/3 to 
2/3 of the SPH smoothing length. Our method has the advan- 
tage that we can readily compute the density in an RT cell from 
the number of particles in the cell and the cell size, whereas 
in the Oxley & Woolfson method the density must be com- 
puted by an SPH sum over the neighbouring particles, and this 
is computationally expensive. 

The method we use to construct RT cells, is implemented 
relatively easily within SPH and exploits information about the 
tree already calculated for SPH purposes. Thus, it does not add 
much to the code running time. It also furnishes us with an 
efficient procedure for finding which RT cell an L-packet is 
in. This is critical for the efficiency of the code, because as 
an L-packet propagates through the computational domain, this 
search procedure is performed many times. 

2.2. L-packet propagation and interaction 

The propagation of L-packets through the computational do- 
main is performed in small steps, 5S , which gradually decrease 
the optical depth r of the L-packet until t = 0, whereupon the 
L-packet reaches an interaction point. This is the most compu- 
tationally expensive part of the code, particularly as there is no 
analytical expression for the density - and hence the opacity 
- along the path of an L-packet. Instead, the density is deter- 
mined by identifying the RT cell through which the L-packet is 
passing. After each step 6S , we use the SPH tree to find which 
RT cell i the L-packet is now in. The search starts from the root 
cell and proceeds to lower level cells, until the required RT cell 
is reached. The search contributes a significant fraction of the 
code running time, and therefore, it is essential to choose the 
step 6S appropriately. We use 

SS — MIN {771 Si, 772 lu T iu ?73 |r|} > (2) 

where 771, 772, 773 are constants in the range (0.1, 1) and deter- 
mine the accuracy; S ; is the linear extent of RT cell i; i, = K\pi 
is the mean-free-path for the L-packet (K\ is the monchromatic 
mass-opacity at the wavelength of the L-packet, and p, is the 
density in cell 7); r is the position of the photon relative to any 
discrete source (i.e. star). 
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Kurosawa & Hiller (2001) adopt a more detailed approach, 
in which they use the path of an L-packet to determine where 
it intercepts the boundaries of the RT cells through which it 
passes. Our method is simpler to implement, it propagates the 
L-packets efficiently (usually the propagation step is on the or- 
der of the local RT cell size), and it is accurate (due to the fact 
that neighbouring RT cells have similar densities). 



3. Tests 

We consider a spherical cloud with uniform density n = 8.6 x 
10 4 cnr 3 , mass M = 100 Mq and radius R = 3.5 x 10 4 AU, 
and represent it with randomly distributed SPH particles. We 
then test the method of constructing the global grid and propa- 
gating L-packets through the grid in three ways, (i) We use the 
thermodynamic equilibrium test, where the cloud is bathed in 
an isotropic blackbody radiation field with a given temperature 
(Stamatellos & Whitworth 2003). (ii) We illuminate the cloud 
with the Black (1994) interstellar radiation field (BISRF) and 
compare our results with the results obtained previously using 
a spherically symmetric grid of concentric cells (Stamatellos & 
Whitworth 2003). (iii) We embed a low-temperature star at the 
centre of the cloud and again compare the results with those 
obtained previously using a spherically symmetric grid of con- 
centric cells (see online appendix). We present the results in the 
following subsections. 



3.1. Cloud bathed in a blackbody radiation field 

We use Af TOTAL = 20,000 SPH particles and set the maxi- 
mum number of particles in an RT cell to N MAX = 60, giving 
^cells = 1,484 RT cells, each containing on average N ~ 15 
SPH particles. The density in each RT cell is plotted as a func- 
tion of radius in Fig. Because the SPH particles are dis- 
tributed randomly, there are statistical variations in the density 
at the ±N~ l/2 ^ +25% level. There are also a number of RT 
cells that appear to be outside the boundary of the cloud and 
to have very low density. In reality, only a small part of each 
of these RT cells is within the cloud and therefore they contain 
only a very small number of particles. These irregular RT cells 
are a side effect of representing a spherical cloud with a grid of 
cubic cells. The same problem was also encountered by Oxley 
&Woolfson(2001). 

10 8 L-packets are then injected from random positions on 
the boundary of the cloud, and with random frequencies and in- 
jection angles, so as to imitate an isotropic blackbody radiation 
field at 15 K (see Stamatellos et al. 2004). The cloud should 
then adopt the temperature of the illuminating radiation field, 
i.e. 15 K. The computed temperature of each RT cell is plot- 
ted on Fig.^3. In the body of the cloud the temperature is very 
close to 15 K, with the error being less than + 3%. The errors 
are higher near and outside the boundary of the cloud, because 
of the presence of the irregular RT cells. These errors can be 
reduced by increasing the number of SPH particles, and/or the 
number of RT cells, and/or the number of L-packets. 
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Fig. 1. A uniform-density spherical cloud, represented by 
20,000 SPH particles, is illuminated by an isotropic blackbody 
radiation field at T = 15 K, using 10 8 L-packets. RT cells 
are constructed with A^ MAX = 60, giving jV cells = 1,484. (a) 
Triangles give the densities of the individual RT cells, the hori- 
zontal dotted line marks the mean density, and the vertical line 
marks the boundary of the cloud, (b) Triangles give the tem- 
peratures of the individual RT cells, the horizontal dotted line 
marks T = 15 K, and the vertical line marks the boundary of the 
cloud, (c) The percentage error in the computed temperature. 



3.2. Cloud embedded in the ISRF 

We use iV TOTAL = 200,000 SPH particles and set the maxi- 
mum number of particles in an RT cell to N MAX = 100, giving 
^cells = 1 1, 862 RT cells, each containing on average N ~ 18 
SPH particles. 10 8 L-packets are then injected from random 
positions on the boundary of the cloud, and with random fre- 
quencies and injection angles, so as to imitate the isotropic 
Black (1994) interstellar radiation field (BISRF). The advan- 
tage of this test is that the BISRF covers a wide wavelength 
range, including regions where scattering dominates (UV and 
optical) and regions where absorption/emission dominates (IR 
and submm). Therefore, we can test the validity of the method 
over a wide range of wavelengths. Fig.|2]compares the result- 
ing temperature profile and SED with accurate results obtained 
using a spherically symmetric grid of thin concentric RT cells. 
The agreement is very good, apart from the irregular RT cells 
at the boundary. 



4. Stars in SPH-RT simulations 

The presence of a normal- or high-temperature star in an SPH 
simulation makes the treatment of radiation transport more 
complicated, firstly because in the immediate vicinity of the 
star the dust is destroyed by sublimation and/or chemical sput- 
tering (e.g. Lenzuni et al. 1995), and secondly because just 
outside this region there are very large temperature gradients 
which cannot be captured by the global grid. 



4.1. Dust destruction radius 

The dust destruction temperature is estimated to be between 
L DEST ~ 800 K and L DEST ~ 2100 K (Lenzuni et al. 1995; Duschl 
et al. 1996) and depends on the assumed dust composition. The 
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Fig. 2. A uniform-density spherical cloud, represented by 
200,000 SPH particles, is illuminated by the BISRF radia- 
tion field, using 10 8 L-packets. RT cells are constructed with 



100 giving JV C1 



11,862. (a) Triangles give the 



temperatures of the individual RT cells and the solid line gives 
the accurate temperature profile, (b) The computed SED (solid 
line) is compared with the accurate SED (dashed line) and the 
illuminating BISRF (dotted line). 

heating rate for an isolated dust grain at distance R from a star 
having radius and surface temperature L+ is 



Q 



nB Y {Ti,) dv , 



(3) 



where r& is the radius of the grain, and Q Y is the absorption 
efficiency. Similarly the dust cooling rate is 



£ = 4nr, 



Jo 



Q v nB y (T)dv, 



(4) 



where T is the dust temperature. Assuming that the dust is in 
thermal equilibrium, i.e. Q — X, and putting Q v oc y (1 i B i 
2), we obtain 



,2/(4+/!) 



(5) 



At the dust destruction radius, T = L DEST , and hence R = R v 
where 



(4+yS)/2 



(6) 



In the optically thin limit, this radius defines the spherical vol- 
ume around the star which should be devoid of dust. If the en- 
velope of dust around the star is optically thick to its own cool- 
ing radiation, it will act to trap the radiation from the star and 
R DEST will be somewhat larger. The values of ^ DEST stipulated by 
Ivezic et al. (1997; see below) take this into account. 

4.2. Steep temperature gradients in the vicinity of a 
star 

The density and temperature gradients close to a star are large, 
and can only be resolved with RT cells whose radial extent is 
comparable with, or smaller than, the dust destruction radius 
R DEST - Unless an impractically large number of SPH particles 
has been used in the simulation, the cubic RT cells of the global 
grid (with size ~ h, where h is the smoothing length of the SPH 




star grid 



temperature 
the star grid 



temperature not 
using the star grid 



Rde; 
Distance from the star 



Fig. 3. (a) The square cells represent the global grid derived 
from the SPH tree, and the concentric circles represent the ad- 
ditional star grid, (b) The dashed horizontal lines show the tem- 
perature computed (inaccurately) using the global grid on its 
own, whilst the solid curve shows the (accurate) temperature 
profile computed using the additional star grid. 



kernel) are much larger than R DEST , as illustrated schematically 
in Fig. 13 Thus the steep temperature gradients in the vicinity 
of the star occupy only a few RT cells of the global grid, and 
they are not properly resolved. Consequently the temperatures 
very close to the star are underestimated. (Further out the tem- 
perature may be overestimated or underestimated, but this is 
less critical.) The upshot is that L-packets which are absorbed 
very close to the star are then re-emitted from a wavelength dis- 
tribution corresponding to an inappropriately low temperature, 
and hence at inappropriately long wavelengths. This means that 
they are able to propagate further away from the star, or even 
to escape from the system altogether. Consequently the SED 
has a deficit of radiation at the short (mainly NIR) wavelengths 
emitted by hot dust, and a corresponding excess of radiation at 
long (mainly FIR) wavelengths. This problem also arises in the 
simulations of Kurosawa et al. (2004). 

To resolve this problem, we construct an extra grid of con- 
centric shells (a star grid) around each star (see Fig. [3}- To con- 
struct a star grid, we must specify: 
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- The dust destruction radius, R DEST , which is determined by 
the temperature T+ and the radius R+ of the star. For sim- 
plicity we put 

, d /AM 3 GM+M+ 
L.=4^=L (-) + , (7) 

where M* is the mass of the star and M* is the accretion 
rate onto the star. The first term in Eq. is the intrinsic 
luminosity of the star, and the second term is the accretion 
luminosity. For a Pre-Main Sequence star we can put R+ st 
3R Q (Stahler 1988). 

- The outer radius R SG , of the star grid, which should be on 
the order of the local SPH resolution, or, equivalently, on 
the order of the size of the local cubic RT cells of the global 
grid. 

- The density profile between R DEST and R sa , for which we 
adopt 

with po fixed by conservation of mass, and p = to 2 . 

Once all the above parameters have been specified, the star 
grid around each star is constructed by dividing up the space 
between ^ DEST and R sc into concentric RT cells of equal loga- 
rithmic width; typically ~ 30 RT cells are sufficient for each 
star grid. The complication that this introduces into an RT cal- 
culation is minimal. There is some overlap between the two sets 
of RT cells (the cubic cells of the global grid and the spherical 
cells of the star grid), but this does not create any major prob- 
lems, as our tests in the next section demonstrate. L-packets 
which are inside the star grid interact with the RT cells of that 
grid, rather than the RT cells of the global grid. 

However, we note that a spherically symmetric star grid 
may not be appropriate when the dust close to a star is con- 
centrated in a disc. This situation can be handled with a simple 
extension to the method we have presented here, and it will be 
discussed in the companion paper (Stamatellos et al. 2005). 

5. Ivezic test 

We test the above method for including stars in radiative trans- 
fer calculations within SPH, using the configuration proposed 
by Ivezic et al. (1997), viz. a star with L* = 2500 K surrounded 
by a spherical envelope with an R 2 density profile. We perform 
the test for envelopes having visual radial optical depths (from 
the centre to the edge of the envelope) ry = 1 and Ty = 100 
(see online appendix). For Ty = 1, we set ^ DEST = 9.11/?*; 
and for Ty = 100, R DEST = 17.677?*, as stipulated by Ivezic et 
al. (1997). For both cases we set R SG = 2QR DEST and R omER = 
10QQR DEST (where R OUTER is the outer radius of the envelope). 
Note that we do not need to specify R+. 

We represent the envelope using 20,000 SPH particles, dis- 
tributed randomly so that they approximate to an R 2 profile. 
(The approximation could be improved by settling the distribu- 
tion to get rid of Poisson fluctuations, but we have not done 
this.) We use the SPH tree to construct a global grid with 



N MAX = 55 (which gives A^ CELLS = 1276 and N = 15). In ad- 
dition, we have the option to introduce a star grid comprising 
30 concentric RT cells, between ^ DEST and R SG . 

The results are compared with the very accurate results 
obtained using a one-dimensional spherically symmetric code 
in Stamatellos & Whitworth (2003), hereafter 'the benchmark 
calculation' . 

We have performed the Ivezic test with Ty = 1 and p - -2, 
using 5xl0 7 L-packets, both without (Fig.0}, and with (Fig. [5}, 
a star grid around the central star. 

Without a star grid around the central star (Fig.|4j the cubic 
RT cells of the global grid in the vicinity of the star have size 
~ 2QR DEST , and therefore they are unable to capture the large 
temperatures and densities there. This has two effects. First, 
the hottest dust is only at ~ 300 K (rather than ~ 800 K in 
the benchmark calculation), and consequently the SED shows 
a deficit of radiation at wavelengths around 8 pm. Second, the 
uniform density in the individual RT cells reduces the net radial 
optical depth, allowing more radiation to escape directly from 

EFFECTIVE 

the central star (specifically t v ~ 0.3); this shows up in 
the SED as an excess of radiation at wavelengths around 1 pm. 

With a star grid around the central star (Fig. qj, the in- 
nermost RT cell (immediately outside ^ DEST ) has radial extent 
~ 0. 1^ DEST , and the high temperatures and densities near the star 
are properly represented. Agreement with the benchmark cal- 
culation is then very good. Small differences are due to statis- 
tical noise in RT cells which do not intercept many L-packets. 

6. Discussion 

In this paper we have presented a method for performing radia- 
tive transfer calculations on individual time-frames from SPH 
simulations. We use the SPH tree to construct a global grid of 
radiative transfer (RT) cells. This grid is constructed so that 
each RT cell contains a number of SPH particles close to the 
number of neighbours that each particle must have in SPH for 
smoothing purposes. The method thereby creates RT cells with 
size on the order of the local SPH smoothing length. 

We test the method using the thermodynamic equilibrium 
test for a uniform-density spherical cloud, bathed in a black- 
body radiation field. The method performs this test well. The 
global grid of RT cells that represent the cloud acquire the same 
temperature as the blackbody radiation field. The temperature 
is slightly underestimated near the boundary of the cloud due 
to the presence of RT cells that contain only a few SPH parti- 
cles (irregular RT cells). The method also performs well for a 
uniform-density cloud illuminated by the interstellar radiation 
field. 

We then examine a system comprising a low-temperature 
star embedded at the centre of a uniform-density spherical 
cloud. The temperature field and the SED are calculated with 
good accuracy. Very close to the star the temperature is slightly 
underestimated due to the fact that the temperature field there 
is not properly resolved. This in turn distorts the SED of the 
system slightly at short wavelengths. The problem is more pro- 
nounced when a normal- or high-temperature star is embedded 
in the cloud. The region very close to the star is then char- 
acterised by a very steep temperature gradient which cannot 
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Fig. 4. Ivezic test with t v = 1, using 5 x 10 7 L-packets. The 
envelope is represented by 20,000 SPH particles and a global 
grid of 7V CELLS = 1,276 RT cells (with N = 15), but no addi- 
tional star grid around the star, (a) The triangles represent the 
temperatures of the individual RT cells plotted against radius, 
whilst the solid line represents the benchmark calculation. The 
global grid on its own fails to capture the temperature profile 
close to the star, (b) The solid line shows the SED computed 
here, and the dotted line shows the benchmark calculation, (c) 
The percentage difference between the flux computed here and 
the benchmark calculation. 





Fig. 5. Same as Fig.|4]but with an additional star grid around 
the star, i.e. a volume devoid of matter inside R DEST and 30 con- 
centric spherical RT cells between R DEST and R sa = 20^ DEST . 
Agreement with the benchmark calculation is now very good. 
Small differences are due to statistical noise. 



normally be captured by the cubic RT cells of the global grid, 
because these RT cells are too large. 

To solve this problem we introduce an additional grid of 
concentric spherical RT cells around the star (the star grid). 
The star grid accounts for the dust-free region near the star, and 
its concentric spherical RT cells are small enough to capture 
steep temperature gradients, so that temperatures are calculated 
correctly close to the star. We test the use of a global grid in 
tandem with a star grid around the star, using the configuration 
proposed by Ivezic et al. (1997), and show that it performs very 
well, even if the optical depth is very large. 

Our radiative transfer code can therefore be applied to arbi- 
trary SPH density fields to produce dust temperature distribu- 
tions, SEDs and isophotal maps at different wavelengths. These 
can then be used to compare the results of hydrodynamic simu- 
lations with observations. In the companion paper (Stamatellos 
et al. 2005) we apply this method to the collapse of turbulent 
molecular cores and the early stages of star formation. 



Acknowledgements. We gratefully acknowledge support from 
the EC Research Training Network "The Formation and 
Evolution of Young Stellar Clusters" (HPRN-CT-2000-00155) 
and from PPARC (PPA/G/O/2002/00497). We would also like 
to thank J. Barnes for using his tree code (available online at 
http://www.ifa.hawaii.edu/~barnes/software.html l. 

References 

Baes, M., Stamatellos, D., Davies, J., Whitworth, A. P, et al. 2005, 

New Ain press, l astro-ph/0504037 l 
Barnes, J. & Hut, P. 1986, Nature, 324, 446 

Benz, W. 1990, Numerical Modelling of Nonlinear Stellar Pulsations 

Problems and Prospects, 269 
Benz, W. 1991, LNP Vol. 373: Late Stages of Stellar 

Evolution. Computational Methods in Astrophysical 

Hydrodynamics, 259 
Black, J. H. 1994, ASP Conf. Ser. 58: The First Symposium on the 

Infrared Cirrus and Diffuse Interstellar Clouds, 355 
Bjorkman, I. E. & Wood, K. 2001, Apl, 554, 615 
Duschl, W. }., Gail, H.-P, & Tscharnuter, W. M. 1996, A&A, 312, 624 
Efstathiou A. & Rowan-Robinson M., 1990, MNRAS, 245, 275 
Efstathiou A. & Rowan-Robinson M., 1991, MNRAS, 252, 528 
Gingold, R. A.& Monaghan, 1. 1., 1977, MNRAS, 181, 375 
Henyey, L. C. & Greenstein, I. L. 1941, Apl, 93, 70 
Hernquist, L. & Katz, N, 1989, ApIS, 70, 419 

Ivezic, Z., Groenewegen, M. A. T, Men'shchikov, A., & Szczerba, R. 

1997, MNRAS, 291, 121 
Kessel-Deynet, O. & Burkert, A. 2000, MNRAS, 315, 713 
Kurosawa, R. & Hillier, D. J. 2001, A&A, 379, 336 
Kurosawa, R., Harries, T. I., Bate, M. R., & Symington, N. H. 2004, 

MNRAS, 351, 1134 
Lefevre I., Bergeat I. & Daniel I. Y. 1982, A&A, 1 14, 341 
Lefevre I., Daniel J. Y. & Bergeat J. 1983 A&A, 121, 51 
Lenzuni, P., Gail, H, & Henning, T. 1995, Apl, 447, 848 
Lucy, L. B. 1977, Al, 82, 1013 
Lucy L. B. 1999, A&A, 344, 282 

Men'shchikov A. B. & Henning Th. 1997, A&A, 318, 879 
Monaghan, 1. 1. 1992, ARA&A, 30, 543 
Oxley, S. & Woolfson, M. M. 2003, MNRAS, 343, 900 
Stahler, S. W. 1988, Apl, 332, 804 

Stamatellos, D. & Whitworth, A. P., Boyd, D.F.A, & Goodwin,S. P. 

2005, A&Ain press 
Stamatellos, D. & Whitworth, A. P. 2003, A&A, 407, 941 
Stamatellos, D., Whitworth, A. P., Andre, P., & Ward-Thompson, D. 

2004, A&A, 420, 1009 
Viau, S., Bastien, P., & Cha, S-H 2004, in preparation 
Whitehouse, S. C. & Bate, M. R. 2004, MNRAS, 353, 1078 
Whitney, B. A., Wood, K, Bjorkman, I. E„ & Wolff, M. I. 2003, Apl, 

591, 1049 

Wolf, S., Henning, T, & Stecklum, B. 1999, A&A, 349, 839 

Wood, K, Lada, C. I., Bjorkman, I. E., Kenyon, S. I., Whitney, B., & 

Wolff, M. I. 2002, Apl, 567, 1 183 
Zucconi, A., Walmsley, C. M., & Galli, D. 2001, A&A, 376, 650 



Dimitris Stamatellos & Anthony P. Whitworth: Radiative transfer in SPH density fields, Online Material p 1 



Online Material 



Dimitris Stamatellos & Anthony P. Whitworth: Radiative transfer in SPH density fields, Online Material p 2 



Appendix A: Supplementary tests 

A.1. Low-temperature star surrounded by a spherical 
envelope 

We use N TOTAL = 200,000 SPH particles and set the maxi- 
mum number of particles in an RT cell to N MAX = 60, giving 
^cells = 13,050 cells, each containing on average N ~ 16 
SPH particles. 10 s L-packets are emitted by a central low- 
temperature 'star' having luminosity L* = 7L and surface 
temperature L* = 30 K. Fig. IA. ll compares the resulting tem- 
perature profile and SED with accurate results obtained using 
a spherically symmetric grid of thin concentric RT cells. The 
global grid used here has difficulty capturing the steep tem- 
perature gradient near the star. As a result the temperatures 
near the star are underestimated and the SED shows a deficit 
at short wavelengths. The temperature is also not very accurate 
near the boundary, due to the irregular RT cells there; and there 
are some small uncertainties in the SED at wavelengths longer 
than 2000 pm, due to the small number of L-packets generated 
at these wavelengths. Otherwise the agreement is good. 




Fig.A.l. A uniform-density spherical cloud, represented by 
200,000 SPH particles, is heated by a central low-temperature 
star with L* = 7L and T* = 30 K, which emits 10 8 L- 
packets. RT cells are constructed with A^ MAX = 60, giving 
^cells = 13,050. (a) Triangles give the temperatures of the 
individual RT cells and the solid line gives the accurate temper- 
ature profile, (b) The computed SED (solid line) is compared 
with the accurate SED (dashed line), (c) The precentage error 
in the SED. 



A.2. Ivezic test with t v = 100 and p - -2 

We have performed the Ivezic test with Ty = 100 and p = 
-2, using 10 7 L-packets, both without (Fig. IA.21 . and with 
(Fig. IA.3t . a star grid around the central star. 

Without a star grid around the central star ("Fig. IA.2L the cu- 
bic RT cells of the global grid in the vicinity of the star are again 
unable to capture the large temperatures and densities there. 
Consequently the hottest dust is only at ~ 200 K (rather than 
~ 800 K in the benchmark calculation). The SED is dominated 
by the emission from this dust, which peaks around 20 pm. In 
addition, the uniform densities in the RT cells reduces the ef- 

EFFECTIVE 

fective radial optical depth (t v ~ 20) and so there is a 



small shoulder at ~ 3 pm, due to radiation which escapes di- 
rectly from the central star. 

With a star grid around the central star (Fig. IA.3t agreement 
with the benchmark calculation is very good. Small differences 
are again due to statistical noise. We emphasize that these cal- 
culations were performed with only 10 7 L-packets. 
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Fig.A.2. Ivezic test with Ty = 100, using 10 7 L-packets. The 
envelope is represented by 20,000 SPH particles and a global 
grid of 7V CELLS = 1, 310 RT cells (with N = 15), but no addi- 
tional star grid around the star, (a) The triangles represent the 
temperatures of the individual RT cells plotted against radius, 
whilst the solid line represents the benchmark calculation. The 
global grid on its own fails to capture the temperature profile 
close to the star, (b) The solid line shows the SED computed 
here, and the dotted line shows the benchmark calculation, (c) 
The percentage difference between the flux computed here and 
the benchmark calculation. 
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Fig.A.3. Same as Fig. IA.2I but with an additional star grid 
around the star, i.e. a volume devoid of matter inside /? DEST 
and 30 concentric spherical RT cells between /? DEST and R SG - 
20/? DEST . Agreement with the benchmark calculation is now 
very good. Small differences are due to statistical noise. 



